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A FULLY DISCRETE VARIATIONAL SCHEME FOR SOLVING NONLINEAR 
FOKKER-PLANCK EQUATIONS IN HIGHER SPAGE DIMENSIONS 


OLIVER JUNGE AND DANIEL MATTHES AND HORST OSBERGER 


Abstract. We introduce a novel spatio-temporal discretization for nonlinear Fokker-Planck 
equations on the multi-dimensional unit cube. This discretization is based on two structural 
properties of these equations: the first is the representation as a gradient flow of an entropy 
functional in the L^-Wasserstein metric, the second is the Lagrangian nature, meaning that so¬ 
lutions can be written as the push forward transformation of the initial density under suitable 
flow maps. The resulting numerical scheme is entropy diminishing and mass conserving. Fur¬ 
ther, the scheme is weakly stable, which allows us to prove convergence under certain regularity 
assumptions. Finally, we present results from numerical experiments in space dimension d= 2. 


1. Introduction 

In this paper, we propose and study spatio-temporal discrete approximations for solutions to the 
following initial boundary value problem for a non-linear Fokker-Planck equation: 


dtu = AP(u) + div(uVU) 

for (t; x) 6 K>o x 17, 

(1) 

n • VP(u) = 0 

for (7; x) e IR>o x dfl, 

(2) 

rt(0; x) = u°{x) > 0 

for a: 6 17. 

( 3 ) 


The domain 11 = (0,1)^^ is the d-dimensional unit cube, P : M>o ->■ K is the pressure function, and 
V : n ^ M is a smooth external potential. We assume that P is continuous, is smooth on K>o, 
and is strictly increasing. Typical choices for the pressure function are P(w) = u (heat equation), 
and P(u) = u"* with some exponent m > 1 (porous medium equation) or 0 < m < 1 (fast diffusion 
equation). To keep technicalities to a minimum, we shall impose some further restrictions on P, V 
and later, see Section 

l.I. Concept of the discretization. The discretization approach presented here has two main 
features: first, it is variational, and second, it is fully Lagrangian. 

Let us start with a discussion of the Lagrangian structure of the evolution: consider the 
Fokker-Planck equation Q as a transport equation for the density u, 

djU = - div (uv[u]), (4) 

with a vector field v[u] : 17 -*■ that depends sensitively on u, namely 

v[u] =-v[d'(u) + V], (5) 
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where h : K>o ^ K is implicitly defined (up to an affine normalization) by P^('r) = rh"{r) for all r > 0. 
In view of the boundary conditions ([^, this representation implies conservation of non-negativity 
and total mass M of the solutions t >->■ Ut = u{t; ■ ), two key properties that are inherited by our 
discretization. 

Moreover, the transport formulation 0 implies that the solution ut at any time t > 0 can be 
written as the push-forward of an a priori fixed reference density u with total mass M on il under 
a suitable diffeomorphism Tt : ->• fl, i.e., 


Ut = = 


u 

det DTi 




( 6 ) 


Canonical choices for u are either the unit density on 17, or the initial condition . The T* are far 
from unique, but there is a distinguished family (T()t>o, namely the one obtained as flow maps for 
the ordinary differential equation x = v[rt(](a;). In other words, one chooses the maps such that 


dtTt = v[Mt]oTj, 


(7) 


which determines the Tj uniquely up to a global change of variables that is given by the initial map 
To- The latter is arbitrary but subject to the consistency condition {^q)^u = 

For the Lagrangian discretization, we approximate the evolution of the maps Tt in a numerically 
tractable form. Our basic idea here is not to restrict the diffeomorphisms Tt themselves to an ansatz 
space, but their time increments. More specifically, given a time step width r > 0, our ansatz for 
the approximation T(^ of the map Tt at time t = nr is the concatenation of n close-to-identity 
diffeomorphisms t^,..., 


rrn _ ,n ,i 


n-1 


■ o t 


-A ’ 


( 8 ) 


The in ([^ are inductively determined inside a fixed finite-dimensional manifold of diffeomor¬ 
phisms. Therefore, each has a simple form, but the complexity of T^ increases rapidly with 

n. 

For our definition of the recursion relation of the incremental diffeomorphisms we build on 
the variational structure of Q. Observe that the evolution Q for the diffeomorphisms T* is 
actually autonomous, since one can substitute Ut = (Tt)^?2: 


9iTt = v[(Tt)#M] oTt. (9) 

Now, it turns out that 0 is a gradient flow on the manifold of diffeomorphisms on 17 with respect 
to the L^(?2)-metric; the respective potential E is given in (16) below. This structure is not a 
coincidence, but is a consequence of the fact that the original equation 0 is a gradient flow of 
an appropriate entropy functional E with respect to the L^-Wasserstein metric; these relations are 
detailed in Section ^2?2\ below. 

The variational discretization of this gradient flow is now performed using minimizing movements, 
which are a variational formulation of the implicit Euler scheme. Given the diffeomorphism T^~^ 
from the previous step, the next diffeomorphism T^ = is obtained by choosing the increment 

as the minimizer of the functional 


1 


— lltoTr^-Tr^ 


+ E(toTA“^ 


) 


( 10 ) 


||2 

over our finite-dimensional set of diffeomorphisms. Thanks to its variational nature, this particular 
discretization approach has — at least — two valuable analytical properties. The first, obvious one 
is the monotonicity of the entropy E resp. E in discrete time. The second, more subtle one is the 
strict convexity of the minimization problem (10), see Proposition for details. This guarantees 
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not only uniqueness of the minimizer, but also efficient numerical solvability. Convexity of (101 is 
also the key ingredient for our stability analysis. 


1.2. Related discretizations from the literature. The general idea to perform a fully La- 
grangian discretization for problems of the type 0-(i dates back at least to the beautiful but 
apparently almost unknown paper by Russo |22| . In fact, there is an earlier short communication 
by MacCamy and Socolovsky m in which a similar idea for use in one space dimension is sketched. 
Later, Budd et al [6] used a variant of that discretization under the name “moving mesh method” 
to track asymptotically self-similar solutions of the porous medium equation — which is 0 with 
P(u) = (m > 1) and V = Q — in one space dimension numerically. None of these schemes has 
been variational, i.e., the gradient flow structure did not play any role. 

The first use of the interrelation between the Lagrangian and the gradient flow character for 
numerical purposes is apparently due to Gosse and Toscani [13] . The method developed there, 
however, was limited to one space dimension, where the Wasserstein metric between probability 
measures on the line is isometrically equivalent to the L^-distance between the corresponding inverse 
distribution functions; the higher-dimensional situation is much more complicated, as discussed 
above. Subsequently, the one-dimensional method has been refined and then applied to various 
further evolution equations, like the Keller-Segel model |S], the Hele-Shaw flow m, the quantum 
drift diffusion equation m, the p-Laplace equation [1] , and even to the isentropic Euler equations 
[25] . For some of those schemes, the discrete-to-continuous limit has been rigorously established by 
the second and the third author [nunj. 

In contrast, the first attempts to combine a fully Lagrangian discretization with a variational time 
integrator for solution of the multi-dimensional problem 0-0 are only quite recent. In fact, we 
are aware of only one fully Lagrangian scheme for solution of non-linear Fokker-Planck equations 
that also respects the variational structure of the underlying gradient flow: this is the work by 
Carlier et al |3|. We emphasize that although their discretization approach has several similarities 
to our own — e.g., it produces a sequence of strictly convex minimization problems — the way they 
exploit the Lagrangian structure is conceptually different from the idea to iterate close-to-identity 
diffeomorphisms from a given ansatz space, as it is done here. It should further be noted that 
some preliminary considerations in the same direction, along with numerical experiments, have 
been made in the thesis of Roessler [^ . 

For completeness, we mention two further works that do not completely ht the context, but are 
still closely related: first, Carrillo and Moll |5] developed a fully Lagrangian scheme for a non-local 
agregation equation; although that scheme is not truely variational, its definition is closely linked 
to the gradient flow structure. Second, in a recent paper by Cavalletti et al |9], the isentropic Euler 
equations have been analyzed in the framework of Wasserstein gradient flows. The construction 
used there can be interpreted as a fully Lagrangian variational solver. 


1.3. Plan of the paper. The motivation and the details of our discretization approach are given 
in Section below. The rest of the paper is divided into two parts: an analytical one, in which 
we perform a rigorous study of well-posedness and consistency, and a numerical one, in which we 
clarify practical aspects of the implementation and report results from experiments. Our results 
in the analytical part are summarized in Section [^ Their respective proofs are given in Sections 
and[^ Section]^ is concerned with the implementation of the numerical scheme. Numerical results 
are then given in Section [^ 


2. Preliminaries 
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2.1. Additional hypotheses and notations. The class of non-negative densities with total mass 

M is denoted by 'P(n). The function q : ^ K is defined by 

= ( 11 ) 

observe that Vq = id, the identity map on For brevity, we introduce 

v(t;T) :=v[(to'I)^M]ot. (12) 

By Diff^, we mean the class of (orientation-preserving) C^-diffeomorphisms T of fl onto itself that 
fix the corners. Note that any T e Diff'*' maps each of the {d- 1)-dimensional boundaries of ft onto 
itself, that is, boundary points are only moved laterally. 

Symmetry assumption: To avoid unnecessary technicalities, we restrict our considerations to 
solutions u with a certain symmetry. Specifically, we call a function / : fl -»• M symmetric-periodic 
if it admits a continuous extension / : ^ M to the whole space which is both even and 2-periodic 

with respect to each coordinate. We further say that / is -symmetric-periodic if additionally 
/ € C^(K‘^). For instance, any function of the form 

f{x) = F{costtxi, ... ,cos7ra;d) 

with a C^-function F : ->■ K is C*^-symmetric-periodic. Smooth symmetric-periodic solutions 

u form an invariant class for the problem in the following sense: assuming that both the 

potential V and the initial datum u® are C'°“-symmetric periodic, it then follows by unique smooth 
solvability of 0-@, see e.g. |23j . in combinaition with symmetry arguments that the corresponding 
solution ut is C°“-symmetric-periodic at any t > 0. 

Variational derivative and gradient: In order to avoid the use of the ambigous symbol SFjSu for 
the variational derivative, we introduce the following convention: for a given functional F, denote 
by DF(m)[(p] the first variation of F at u in the direction of (the smooth function) (p. Assuming 
that DF(u) extends to a linear and continuous functional on L^{9), where 9 is some given measure 
on ft, let JL 2 (e) [DF(u)] e L^{9) be its Riesz-dual, i.e., 

Ax) [D^(m)]} (a;) d6»(a:) = DF(m) [p] 

holds for all smooth p. 

2.2. Analytical background — two gradient flow structures, Lagrangian representation. 

To motivate the variational discretization approach, we briefly discuss the gradient flow structure 
behind which is actually two-fold. 

2.2.1. First gradient flow structure. One possibility to write ([l|)-([3|) in variational form has been 
worked out in Otto’s celebrated paper [201 • Observe that v in^^ has the form 

v[m] = -VJL2(da;) [DE(u)] , (13) 

where E is the following relative entropy functional 

E(u) = f h(u{x))dx+ f u{x)V{x)dx. (14) 

Jn Jn 

This makes solutions t^Ut curves of steepest descent in the potential landscape of E on the space 
V{Vl) of densities with respect to the F^-Wasserstein metric >V 2 . We refer the reader to for 

an introduction to Wasserstein metrics and associated gradient flows. 

The main analytical insight from this gradient flow representation is that E is a A-convex func¬ 
tional in the metric >V 2 m- More precisely, it is A-convex along generalized geodesics in the sense 
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of [2]- Here A € M coincides with the modulus of convexity of V, i.e., > Al. This makes 

qualitative features of the non-linear Fokker-Planck equation Q accessible to the abstract theory 
of A-contractive gradient flows, with important consequences on the large-time behavior |20j . 


2.2.2. Second gradient flow structure. An alternative variational formulation directly for the trans¬ 
port representation Q&([^ was given by Evans et al in [T^], see also 011]. Indeed, it turns out 
that Q is again a gradient flow, namely on the space of diffeormorphisms on H with respect to the 
Hilbertian structure induced by L^{u). That is, 

v[T#u]o‘I=-J^.(^)[DE(T)], (15) 

for the potential E (on diffeomorphisms T) that is induced by the relative entropy E above: 


E(T) = E(T#u) = ^ 


^ det DT j 


+ 1 / 0-1 


du. 


(16) 


^>0 


above is obtained from h via h*(s) = sh(l/s). Equality of the vector 


The function /i* 

fields V in (13) and in (|15|) is verified by a short calculation, see Lemma 18 


Remark 1. In [12j . the following more explicit representation o/v[T^m] has been derived: 




(Vl/)oT. 


(17) 


Here the divergence div acts on the entries in each row of the matrix (yielding a column vector), 
and denotes the cofactor matrix of A, that is (det A)(A“^)^ = A^. The calculations leading to 
(|17|) are given in the proof of Lemma 19 


We remark that, despite the fact that E is A-convex with respect to the L^-Wasserstein metric 
W 2 , the functional E is poly-convex, but has no other useful convexity properties with respect to 
L\u). 


2.2.3. Comparison of the two structures. Although both gradient flows represent the same Fokker- 
Planck equation, they have different analytical features. The first (“Eulerian”) approach is set 
on the nice space of positive densities, with a A-convex functional, but needs the complicated 
Wasserstein metric. The second (“Lagrangian”) approach uses the convenient L^-metric, but on 
the more complicated space of diffeomorphisms, and lies outside of the setting of contractive flows. 
For our discretization purposes, we rely on the second formulation, but we modify the variational 
problem “in the direction of the Eulerian approach” in order to profit from the convexity there. 
This is made precise in Proposition]^ 

2.2.4. Lagrangian representation of the solution. The following Lemma shows how to construct an 
adapted family of push-forward maps for a Lagrangian representation of a solution to the PDE 
problem 0-il)- It is essentially an adaptation of the result from m Section 3.2] to our situation. 

Lemma 2. Let u : [0, t] x H -> IR>o be a smooth symmetric-periodic classical solution to (j^ . Define 
t: [0, t] X H ->• as the unique solution to the initial value problem 

dAcr = v[ucr] o to-, to = id. (18) 

Then and consequently, i(.) is a solution to daia = v(tcr;T), where T e Diff^ is 

arbitrary with = uq. 
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Proof. Let ip : ■ 

we have that 


be a -symmetric-periodic test function. On the one hand, since 
daUa = ^P{Ua) + V • (UcrW) = -V • {UaV[Ua]), 


dcr 


J' (pUadx= J (pdaUadx=-J' ' v[Mcr] Mcr dx. 


On the other hand, recalling the definition of O in ([T^, 


dcr 


J (f {ta)#Uodx = -^ J {(f o icr) Uo dx = - J (o ter) ' ( v[m<,.] o ) Mg dx 
= -J Vyj- v[m<^] dx. 


Subtracting those equations from each other yields 
d 


dcr 


J (p[ua-{ta)#Uo]dx = - J v[u^] • [u^r “ (tcr)#Uo] dx. 


Choosing (p = 2[u^ - (tcr)^Uo] yields further that 
d 


dcr 


J [ua-{ia)#uof dx = - J v(u„) • V- (t^)#Uo]^ dx 

= J V-v(Ucr)[uo--(to')#Uo]^dx< A J [Uo- - (tcr)#Uo]^dx, 

where A is a uniform bound on the smooth function 

V • v[Ucr] = Ah'{Ua) + AF. 

Since (to)#ito = id# ug = mq, Gronwall’s lemma yields that (tcr)#uo = for all cr e [0, r]. 

3. Discretization 


□ 


3.1. Discretization in time. Below, we motivate our choice of the time discretization by dis¬ 
cussing a semi-discretization of 0-©- Specifically, consider approximations and 

for a given time step r > 0 by means of the minimizing movement scheme: initially, choose e V{Pl) 
as approximation of and T° such that (T°)#?x = vff. Then, define inductively 


= arg mm 

lie'P(Sl) 


-W2{u,K~^Y+¥.{ u) 

It 


A" = arg min 
'X6Diff+ 


1 

2t 




(19) 


By some formal considerations — which have been made rigorous under certain technical hypotheses 
in [ 5 ] — one concludes that the equivalence between the gradient flows is preserved under this 
time-discretization. More precisely, each time step in (19) is moderated by a close-to-identity 
diffeomorphism 1” of fl which is such that, simultaneously, 

n /rn\ n-1 __j n-n pn 

= (4 and = 4 o 

hold. This i" is determined as solution to the corresponding Euler-Lagrange equation 

t= id+rv(t;T”~^), 


( 20 ) 


with the abbreviation v(t;T) defined in (12). One concludes that in each time step 

n, and also that = ||T^ In fact, is an optimal transport maj[^from 

^That is, is a solution to Mongers problem of finding a diffeomorphism t of with such that 

Jq |t(a:) - x\'^u‘!^~^{x) dir is minimal. 
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to m". Thanks to the general theory [5], optimal transport maps t: ^ 17 are those which can be 

written as the gradient of a convex function, i.e., i = for a suitable convex (/): 17 ->• IR. 


Remark 3. We give a heuristic argument why solutions to the Euler-Lagrange equation ( 20 ) are 
optimal transport maps, at least for sufficiently small time steps t. This is not completely obvious 
sinee on the right-hand side of (20), we have the composition of the gradient vector field v with i 
itself; in general, such a composition is not necessarily a gradient again. Fortunately, ( 201 can be 
equivalently written in the form 


i ^ = id-Tv[(to T" 


Thus, if t is a solution, then its inverse is the gradient of some potential / ; 17 ->• K. If this f 
happens to be convex (which is a reasonable guess for sufficiently small t > 0), then it follows by 
elementary convex analysis that t itself is the gradient of a convex function as well, namely of the 
Legendre transformed f*. 

For example, if v[u] = -Vl^ with > A > -oo, then (20) has the unique solution t = 

(id+rVl^)”^. Provided that 1 + tA > 0, this i is the gradient of the Legendre transform of the 
convex function f{x) = \\x(f +tV(x). 


One can therefore restrict the second minimization in (19) — a priori — to maps of the form 
T = toT”“^, where t is an optimal transport. Thus, the second minimization problem in (19) reduces 
to 




■n-1 


with 


t" = argmin 




id 


Il2((T- 


)#“) 




( 21 ) 


where, by abuse of notation, we wrote i 6 VX to indicate that the minimization runs over all maps 
t; 17 ->• 17 can be written in the form t= V(f> for a suitable cj) from the convex affine set 

X:= ; (j) is convex, V(()(17) £ 17}. 

For further reference, we introduce the “regular” subset 


€ X n (7^(17); inf det S/^cj) > 0, ^ - q is symmetric-periodic|, 

recalling the definition of q in I©- Notice that a consequence of </)- q being symmetric-periodic is 
that 


n ■ (Vt/i - id) s 0 on (917, 


which in turn implies that Vcj) maps onto 17, and that points on the four boundary segments only 
move laterally (and the corners stay fixed) under Vcf. In the same way, we understand t e VX^; 
notice that such a t is an optimal transport map and lies in Diff'^. 

The reformulation (21) does not only reduce the complexity of the variational problem (19) — 
by moving from diffeomorphisms to scalar functions — but restores A-convexity: 


Proposition 4. For each T e Diff^, the functional i (->• E(to T) is X-convex on VX in the usual 
(flat) sense, with A being V’s modulus of convexity. 


The proof is given in Section [O] 
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3.2. Discretization in space. Spatial discretization is achieved by restriction of the potentials 4> 
for the close-to-identity transformations t” to a finite-dimensional convex affine set X 5 c X of ansatz 
functions, where S is an indicator for the spatial resolution. Specifically, we assume that a system 
{&k}k6i of C°°-symmetric-periodic ansatz functions bk : -> K is given, that forms an orthogonal 

basis of L^{dx). To each ^ > 0, a subset c X of indices k is associated such that Is' £ whenever 
i5' > (5 > 0. Now introduce 


TX^ = \ip= Zk&k ; ^k 6 K [, 

[ k«sla J 

the set of linear combinations of elements in {Vbklkeij- From that, we define 

Xs = {(j) = c{ + tp(j) € TX5} n X, 

and accordingly XJ. 

Remark 5. A very natural choice for the bk are Fourier modes, so that TX^ is a space of Fourier 
polynomials. This is what we actually use for our numerical experiments, cf. Section^for details. 

For simplicity, we combine the time step r > 0 and the indicator 6 into a single discretization 


parameter A = (t;(5). Now (211 is replaced by 


rrn_Ln rri 
Aa - Ca ° a. 


A J 


with = arg min 
RVXt 




( 22 ) 


Finally, in accordance with ([^, the approximate Lagrangian map at time step n 

is now obtained as concatenation of the n optimal transports € VXJ. The associated 

“discretized” densities are then given by 


“A = 


n=0,l,2,. 


= {XY)#u. 


(23) 


4. Summary of analytical results 
The fully discrete scheme is well-posed: 

Theorem 6. Assume that s s‘^^^”^/i*(s) is not integrable at s = Q, or, equivalently, that r i->- 
is not integrable at r = + 00 . Let the time step t > 0 be sufficiently small to verify 

X + T-^> 0. 


Then the minimization problems in (22) are convex in each step and can he uniquely solved 


inductively. The nth iterative map belongs to XJ and is the unique solution to the following 
system of Euler-Lagrange equations. 


(i^_v(t;Xr'),VV') 


= 0 




for all ip € TX^. 


(24) 


The proof is given in Section [0| 

The next result is concerned with the continuous limit r ->• 0 and i5 -> 0. For a special choice of 
Fourier ansatz functions bk — cf. Section|^for details — we have the following convergence result. 

Theorem 7. Let w- [0,r] x 12 -> K>g be a classical symmetric-periodic solution to ([^-([^. Choose 
e Diff^ such that (X^),^?! = , and define the discrete approximation ua by means of the 

scheme (21), with (23). We make the following stability assumption: the densities are are 
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A-independently bounded in and positively bounded from below, that is, there are constants 

B, P > 0 such that 


IIuaIIc'^ < -B, > P for all T > 0 and n € N with nr < T. 

Then the discrete solutions ua converge to the smooth solution u in Wasserstein, 

yV2(UnT,UA) = 0(T-t- ps), 

where p(.) : (0,1) ^ M>o is an increasing function with ps ^ 0 as S ^ 0, see Lemma 
The proof is given in Section [673| 


11 


(25) 


(26) 

for details. 


Remark 8. • Note that the hypotheses and conclusions are formulated in term of densities, 

not in terms of the transport maps by which they are induced. With additional technical 
effort, an approximation of the transport should be possible, using the techniques from . 

• Since we require (25) in particular at the initial time, it follows from the regularity theory 
for filtration equations that u is -symmetric-periodic. 

• The assumption of A-uniform -bounds on the discrete approximations is certainly a 
strong hypotheses to conclude convergence in the Wasserstein metric. Unfortunately, our 
only stability result is the one related to the convexity of the minimization problem {22), 
which provides an estimate in the Wasserstein metric, see Lemma [TSl but no better. A 
much higher degree of regularity of the densities u\ is needed to control the approximation 
error for the velocity field by means of Lemma \lf\ and Lemma \11\ 

• In any case, by standard interpolation arguments, one easily deduces approximation in better 
spaces, like 


IWriT - Ua\\c‘>‘ 0 

as T ->■ 0 and J ->• 0 simultaneously. The estimate on the approximation error above, 
however, will not be linear anymore in r and ps. 


5. Proof of well-posedness (Theorem 

5.1. Generalized convexity. As a preliminary result of somewhat independent interest, we show 
the convexity property stated in Proposition!^ First, we give the formal argument, which is short 
and instructive, following [191 Theorem 2.2]; afterwards, we sketch the rigorous argument. Here 
“formal” just means that we consider t E(to T) as a functional on the regular space The 

extension of the argument to all of VX requires additional technical effort. 


Proof of Proposition^ for regular arguments. Fix T € Diff^. Let toAi ^ be given. Then ig := 
(1 - s)io + sir belongs to VX^ as well, for every s € [0,1], Therefore, detDL is a continuous and 
strictly positive function on H, and so the following calculations are justified: 


JQ \ u ) Jn 


(detDls) o T 
u/ det DT 


F o (Is o T) du 

f (F o tg) o ddu 

Jn 


X'‘*( 
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Since io = V(/>o and ti = V^i for convex functions 6 it follows that is = V<ps for the convex 

function = (1 - s)(/>o + s(j)i 6 We have 

det Dts = det ((1 - s)V^0o + sV^^i) = det V^</!)o det ((1 - s)l + sM), 

with the matrix M = (which is symmetric and positive definite, and hence diagonal- 
izable with only positive eigenvalues, at every x e il. The eigenvalues of the convex combination 
(1 - s)l + sM are the respective convex combinations of the eigenvalues of M with one. Since the 
determinant is the product of the eigenvalues, it follows — see, for instance, m for a proof — that 

s 1 -^ ( det Dis) ^ 

is a concave map. Since z i->- h^{az) is a convex and decreasing map for each a > 0, it follows that 
s 1 -^ /i*( detDts(x)/T:^i2(a;)) is a convex map as well, for each x 6 O, proving convexity of the first 
integral above. Concerning the second integral, it suffices to observe that, for each fixed x € fl, the 
map s y o ts(x) inherits the A-convexity of V. □ 


In the general case with t € VX, a rigorous argument can be obtained along the lines of [5J 
Proposition 9.3.9]. There, it is proven that for any given optimal transport maps to,ti e VX, the 
functional E is A-convex along the generalized geodesic s Us := ((1 - s)to + sti)^{‘I^u) in the 
Wasserstein space. By definition — see [H Definition 9.2.2] — this means that the scalar function 
s <-* E(us) is A-convex on [0,1]. Since 


E{us) = E(((1 - s)io + sti) o T) 

by definition (16) of E, the proposed convexity statement follows. 
Corollary 9. For any T € Diff^ and t°, € VX^, 


-(ti - t°, v(ii;T) - > AUfi - (27) 

Proof. By Proposition above, the functional t E(t o T) is A-convex on VX with respect to 
and by Lemma [l^ its gradient in that space is given by -v(t;T). Therefore, a Taylor 
expansion yields that 

E(ii o T) > E(t° o T) - (fi - i°, v(t0; ^ ||t' - 


The same is true with the roles of and exchanged. 


Addition of both inequalities provides 

□ 


5.2. Well-posedness of the time iteration. We are now in the position to prove that the 
minimization problems in (221 can be solve iteratively w.r.t. n = 1 , 2 ,.... 


Proof of Theorem^ Fix some T€ Diff^. We consider the functional $ : VX^ ^ Ku {+oo} given by 

$(t) = ^|lt-id|li.(^^^)+E(toT). 

Clearly, $(t) is finite for all t € VXJ, but we cannot rule out $(i) = +oo for general i € VX^. 

Now, since X^ is a convex subset of X, also VXa is a convex subset of VX. Clearly, the finite¬ 
dimensional set Xi is closed and bounded, and therefore compact. Thanks to Proposition]^ the 
restriction to VX^ of the map 1 1 -^ E(to T) is A-convex. Clearly, the map 1 1 -^ ||lt- id is 

1-convex on VX 5 , and therefore, $ is a (A + r“^)-convex functional. In particular, <I> is uniformly 
convex for 0 < r < (-A)"^ if A < 0, and for arbitrary r > 0 if A > 0, respectively, independently of T. 
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Observe further that $ is bounded from below, since h* is non-negative and V is smooth. 
Consequently, $ possesses a unique minimizer i* € VX^. It remains to verify that t* € VXJ. 

Thanks to the properties of the basis functions bk, we have that t* e C^(0), and that t* - id 
is symmetric-periodic. It only remains to verify that Dt* is uniformly positive definite on O. For 
that, it suffices to show that f : fi ^ R with f{x) = detDt*(a;) is a strictly positive function. By 
the properties above, / e (7^(17) is symmetric-periodic. Further, since i* e VX is the gradient of 
a convex function, we immediately have that / is non-negative. Towards a contradiction, assume 
that a positive lower bound does not exist. By continuity, there must exist a point a; e 17 such that 
f{x) = 0. Next, non-negativity of / implies that V f{x) = 0; this is true even at the boundary and 
in the corners of 17, thanks to / being symmetric-periodic. From here, C^-smoothness of / implies 
that 


f{x)<c\x-xf for all a; e ]Bj.(a:) n 17, 

for suitable constants c > 0 and r > 0. Let us show that this implies <I>(t*) = +oo, contradicting 
minimality of t*. As a C^-function, / is bounded on 17, hence a := = inf(it//) > 0. Since 

furthermore h* is a positive and decreasing function, we have the following estimate: 

- aax 

a ) 

j/"Mp = ^ h*(c77/a)77‘^/^“M?7. 

By hypothesis on the non-integrability of s at the origin, the last integral is infinite. 

□ 



6. Proof of weak convergence (Theorem 

6.1. Basis functions and their properties. For the proof of Theorem we assume that the 
family of basis functions = {bk}ksi is given by a tensor product of of Fourier modes. More 
precisely, we take the index set I = Nq \ {0} and define the ansatz function bk for the index 
k = (fci,... ,fcd) by 

i2df2 

&k(a;) = cos(fci7ra:i)-"COs(fcd7ra:d)- (28) 

7r|k|2 

For the finite index sets Is, we take 

= {k € I: |k|oo < K} with AT 6 N such that KS € [1; 2). 

Lemma 10. The family 55 = {bkjkei is a complete orthogonal system in L^{dx). The derived 
family V55 = {Vbkjkei is an orthonormal system in L^{dx). 

Proof. The first claim is a classical result. It can be proven using the two facts that 

• every function / e L^(da;) can be approximated (in norm) by sums of functions of the form 
g{xi,X 2 ,...,Xd) = gi{xi)g 2 {x 2 )-"gd{xd), and 

• the cosines-modes cos{rmTx) for m = 1,2,... form a complete orthogonal system in L^((0,1)). 
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To prove the second part, fix indices k,£ el. Since bk and b^ are symmetric-periodic, the following 
integration by parts is justified: 


(Vbk, Vb^)i2(d3;) = Vbk(a:) • Vbf(ai) dx = - bk(ai) Ab^(x) da: 

= cos{kiTTXi) cos(£i7ra:i) dxij 


= <5k,f 


= 7 r ^|£|2 f bk(a:)b^(a:) da; 

Jn 

cos{kdTTXd) cos{£dTTXd) da;^ 


For the last step, we have used the orthonormality of the cos-functions on the interval [0,7r]. □ 


Thanks to the orthonormality, the L^(da;)-orthogonal projection II^ of vector fields v : ^ 

onto the span of VSs has the following easy representation; 

= XI '^^k)L2(da,)Vbk. 

kslj 


Lemma 11. There exists an increasing function /9(.) : (0,1) ^ IR>o with ps ^ 0 as S ^ 0 such that 
the following is true. Suppose that v 'Tl -* is a vector field, which admits a decomposition in the 
form V = V'0 + C) where '0 : ->• K zs a -symmetric-periodic function, and f is a -symmetric- 

periodic remainder. Then 

Ik-n5[ii]|lc72 < Wyf^Wc^ps + 24 ||c||c 3 - (29) 

Remark 12. The proof provides a rough upper bound on p in the form ps oc 


Proof. Thanks to the linearity of If^, we have that 

IIu-Hi[z;] 11^2 < ||Vz/’-n5[Vz/;]||c2 + ||C - [C]||c2 , 

and thus can estimate the terms related to Vip and to f separately. For C, we simply use the 
i5-uniform continuity of II^ as a map from to C'^(n), 

lie- n5[C]||c2 < IICIIc2 + l|n5[C]|lc2 ^ IlClIca + IlClIca ■ 

For the other estimate, we may assume without loss of generality that if is of vanishing mean. 
Define the Fourier representation ips of tp with accuracy 6 by 

i’s = Y. (30) 

kels 

One has that Vz/'a = n 5 [Vz/;] since 

S/Ips= XI (V', l|k||2bk)L2(da;)Vbk = XI (V', “ Abk)L2 (d,,) V bk = XI (^^>'^^k)L2(d;:c)Vbk = n 5 [Vz/)]. 

keXg keX^ keXg 

Therefore, the claim clearly follows if we can show that 

\\pj-tps\\c3 <Cps\\ip\\ci. (31) 


Clearly, the map tp ips is linear. It is easily seen that any C^-symmetric-periodic function 
Ip can be approximated in (7^(0) by sums of functions g that factorize in the coordinates, i.e., 
g{xi,X 2 , ■.. ,Xd) = gi{xi)g 2 {x 2 )---gd{xd), with each gj € C'"^(IR) being even and 2-periodic. Hence it 
suffices to prove (31) for such product functions g in place of ip. Performing the projection (30) on 
g yields again a factorizing function: 


gsixi , ..., Xd ) = gi , s { xi )--- gd , s { xd ), 
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where the individual terms in the factorization are given in the classical way, 

K 

9j,s{y) = / ^xiu, z)gj{z) dz with the kernel rx{y,z) = Y^cos{kTTy)cos{kTrz). 

Note that our particular choice of the index set Is is important here. Comparing corresponding 
partial derivatives on both side of (31), it is now easily seen that the proof of that estimate reduces 
to showing the following relation in one space dimension; 


sup 

y 


dlTK{y,z)f{z)dz 


< CK , 


(32) 


for £= 1,2,3 and all / € C'^(IR) that are 2-periodic and even, with a constant ck that tends to zero 
as K ^ oo. This is covered by standard estimates on the approximation of smooth functions by 
Fourier series. For instance, [TSl Theorem 2.12] yields ( |3^ with ck K . □ 

Lemma 13. Letu ■ [0,r]xr2 ->• K>o, = u{a] ■), be a smooth symmetric-periodic classical solution 
to 0- ( pl w ith uq = u\ ^ o.nd t the flow map from Lemma^ Under the hypotheses of 

Theorem there is a constant C > 0, independent of r and n, such that 

It -n^pr] 


< C{t + ps). 


(33) 


C2 


Proof. Thanks to the regularity hypothesis (25), assuming that r has been chosen sufficiently small, 
one has the bounds 

||Ma||c5<2i?, tier > /3/2. 

Hence the associated cr-dependent vector field 

Va = v[i2£,] = -V(/l'({tcr) + V) 

is bounded in C''^(H), uniformly in cr e [0 ,t]. Consequently, also its flow map t, given by 


(34) 


dArr = Vn- o L 


to = id. 


(35) 


is (T-uniformly bounded in Lemmaguarantees that (iT)#u\ ^ = Ut, and that 

o L = v(C;T”“^). (36) 


Writing out (18) in integral form yields 

tr = id+ r 

Jo 


Ver o dcT = id +TV,- o + r^Ci 


(37) 


with the remainder term 


Ci = - 

T Jo 


T Jo 

which can be estimated as follows: 


j - (VaOf^-W°tT)dCT 

/ f ° ta'dcr'dcr 

T Jo J(7 da' 

r Jo J ) o J' dcr' da, 


||Ci||c3<^ sup [(||v,,||c4||v<,||c3 + ||v^||c3)(l + IILIJ 3 )] < 5' 

0<f7<r 


( 38 ) 
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where A is a combinatorial factor that originates from the repeated application of the chain rule, 
and B' can — in principle — be calculated from B and j3. Moreover, we have that 


with remainder 


that is controlled in the form 


o +TC 2 , 

C2 = ^ (DVr'VCT)oio 


(39) 


dcr. 


||C2||C3^^ sup [||Vr||c^||v,^||c3(l + lltalica)] ^ -S' 

0<(7<r 


(40) 


We conclude from (37)&(39) that 


t,- - id 


= vi/> + C, 


where ip = -[h'{ua) + V], see (34), and ( = t{(i + C 2 ) is controlled by means of (38)&(40). With 
the help of Lemma [TTj we thus obtain the desired estimate (33). □ 

6.2. Weak stability. 

Lemma 14. Assume that e and write u = e (7^(17) with some T e Diff'*'. Assume 

further that A1 < DF < A“^l for j = 1,2, and X<u< with some A > 0 uniformly on 17. Then 

there exists a constant L — expressible in terms of X, OLfid the norms I|t^llc2, ||f^||c2 and 

111111(73 alone — such that 

||v(ti; T) - v(t2; T) Wmu) < L\\l^ - t'llc^ • (41) 

Proof. Introduce the linear interpolation ig = si^ + (1 - s)i^. Clearly, s i->- C is a smooth curve in 
(7^(17), and we have for all s e [0,1]: 

||ts ||(72 < max(||t^||( 72 , ||t^||( 72 ), A1 < Dtg < A ^1. (42) 

Hence the curve di: [0,1] ^ L^{u) with 

4'(s) = v(ts;(r) = -v[/i'((ts)#?l) + C] o tg = -(Dtg)~^Vh' j - VF o tg 
is continuously differentiable, with derivative 

4-'(s) = (Dtg)-^(Dti-Dt2)^(Dtg)-^Vh'(^^) 

+ vVotg-(i^-t^). 


Above, the tg are differentiated at most twice in space. Using (42) and the fact that r h'{r) and 
r rh"{r) are smooth with bounded derivatives for r € it is easy to conclude that 

'I''(s) : 17 ->■ is actually a continuous function for each s e [0,1], with 


|4/'(s)|<L||ti-t2||c3. 


Integrate this inequality with respect to s 6 [0,1] to obtain (41). 


□ 
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Lemma 15. Assume that V ,V 6 VXJ satisfy the set of estimates 


F - id 


- v(t‘ 








(43) 


for all 6 TX^, with constants pi,P 2 ^ 0. Then 


l|t^ - t^llL2(T#a) ^ ^ ^ 


(44) 


Proof. Subtracting the respective inequalities (|43| for the same fj € TXa from each other yields 

^ {p^ + P^)ll Vl/'||L2(<J^fl). 


tl-t2 


-(v(t^T)-v(t^T)),ViA| 




We choose tp e Xs such that t^ - = Vi/i; with Corollary of E’s A-convexity, we arrive at 


{p^ + P^)\^ - > -||f^ - - (t^ - t^, v(t^;T) - v(t^;T)} 








Division by the norm of yields (44). 

6.3. Weak convergence. 

Proof of Theorem^ We shall prove the following estimate, 

y^2{UnT,u'ff) < (l + C>(T))W2('U(„-l)r,MA"^) +TC>(r + P5). 


□ 


(45) 


The result (26) then trivially follows by induction on n. 

Hence fix some n e N. Recall that where he have introduced u := for brevity. 

Define further the flow maps t(.) as in ( [^ , and the density u* = HapT-j^u. Now we apply the 
triangle inequality as follows: 


m{Unr.u'l)<m{Ur^r.hr)+m{hr,K)+y^2{K.y‘l), 


(46) 


and we shall estimate the three terms on the right-hand side separately. The Hrst one is easy: since 
([^ defines a A-contractive gradient flow in W 2 , 

>V2(u„,,i2x) < e-^">V2(M(„_i),,«l-i). 


To estimate the second term on the right-hand side of ( [46| , recall that Ut = (C)^!! and ft* = 
n 5 [C]#M. We thus have, with the help of (33): 


W2(ft^,ft*) < lie-n5[L]||i2(*) < lie - n5[c]||co <t 


ir - n5[L] 


< Ct{t + ps). 


C2 


Finally, to control the third distance in (46), we use that ft* = n, 5 [tT-]^'u and that u\ = so 

that 


W2(ft:,0< ||n5[L]-tilU2(,i). 
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Here we wish to apply Lemma 15 above, with = n 5 [tT] and Clearly, p 2 = 0. To estimate 

pi, we provide a bound on each of the three terms on the right-hand side of 


n^pr] - id 


-v(n,rc];Tl-i) 


L2(*) 


r II tr||^ 2 (ii) 

+ ||v(C;T-i)-v(n5[C];T^-i)j|^,^.^ 


For the first term above, we have, again by (33), that 


1 


- ||n5[tT] - ir 


L^u) 


< C{t + ps). 


To estimate the second term, recall that the relation (36) holds, and therefore (37) yields 
ir - id 




L^{u) 


^ IkClU^cfi) ^ 'tIICIIco ^ 'r||Cllc3 ^ b ' t . 


In the final step, we use the Lipschitz continuity (41) of the map v( • from to L^{u), and 

apply (33) again: 

||v(C;rA-i)-v(n5[C];(ri-i)||^,^.^<L|rC-n5[C]||c2<LCr(r + p5). 

This concludes the proof of (45), from which the final result (p^ can be easily deduced. 


□ 


7. Implementation 

In practice, performing one time step in the iterative scheme means to find the unique minimizer 
€ VXJ of the convex functional 


Jn 


2 t Jq 


uAx + 


/ hj -^^(detDt)oT”-i +I^(ioT^-i)uda;. (47) 

»/ o \ u j 

Although this minimization problem is finite-dimensional, its numerical solution is not immediate, 
since the integrals cannot be evaluated explicitly. Even if u and /i* happen to be “nice”, the 
main difficulty remains, namely the appearance of which is an (n- l)-fold concatenation of 

transport maps from VXJ. 

Therefore, a further simplification of the minimization problem is needed. Our approach is to 


replace the integrals in (47) by suitable quadrature rules. That is, the integrands are evaluated 


only at a hxed number K of quadrature points Xk e H, and these values are then summed against 
given weights uJk > 0. These are chosen as follows: Fix two positive integers Ki and K 2 . Decompose 
n = [0,1]“^ into Kf-mduiij cubes of equal size. It is sensible to choose Ki proportional to the number 
£ of modes per direction. On each of these cubes, introduce Ar|-many quadrature points using a 
one-dimensional quadrature rule on a tensor product grid. 


The “complicated” function ^ appears in two ways when using quadrature on (47): first, in 
the positions := and second, in the determinants := det Fortunately, 

these quantities can easily be determined by iteration, using that = t a ° 

^%,k = ^k, ^A.fc = iA(^A.fc)i 

^A,k = ^A,k = det BeUxXl) ■ 
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Notice that this iteration is explicit in the sense that the quantities x\ f. and j, can be computed 
(in order to determine in the next time step) after the minimizer ^ in (47) has been obtained. 
The functional in (47) now attains the form 

:^n-l \ 

detDt(xrfe) I +y(t(x2rfe)) u{xk)u}k, (48) 


1 


u u 




l{Xk) 


with the quantities I and that have been calculated in the previous time step. Observe 


that in (48), the reference densities u{xk) are constants, and that and detDt(x^^) are 

easy to compute for any given t € VXJ. 


Lemma 16. The minimization problem (48) is strictly convex inte . There exists a minimizer, 
which is unique if it belongs to VXJ. 

Remark 17. In contrast to Theorem it is not a conclusion that the minimizer t* belongs to 
VXJ. The problem is that Dt* might in principle degenerate at any point x € SI except for the Xk- 

Proof. Convexity follows along the same lines as in the proof of Proposition Indeed, the core of 
the argument is the verihcation of the pointwise convexity of the integrand in the definition E. In 


(48), the integral is replaced by a summation. To obtain the uniqueness, it suffices to observe that 
XJ is in the interior of the convex set XJ. □ 


8. Numerical results 


We present the results of a numerical experiment in space dimension d = 2 and for a quadratic 
pressure function P(s) = s^. Even though this is a borderline case for which Theorem does not 
guarantee the well-posedness of the time-discrete iteration, our scheme works well and produces 
positive bounded fully discrete solutions. For the quadrature rule in (48), we choose Ki = 2K and 
K 2 = 2. Thus, n = [0,1]^ is decomposed into 477^ squares, with four quadrature points Xk per 
square. The nodes Xk and the weights ujk are chosen according to the Gaufi quadrature rule. 


8.1. Reference solution. In order to numerically estimate the order of convergence of our scheme, 
we use a reference solution Mj-ef that is produced by a highly resolved FEM method with fully implicit 
time stepping. As ansatz functions (fk-, we choose tensor products of hat functions with respect to 
a standard square lattice. The nth time iterate u"g£ of the reference solution is then obtained by 
solving 

- — —^^(pkdx= f +u'\/V]-'\Tpkdx for all A: (49) 

T Jn 

for u=Y,k UkTk- We use 400 lattice points in both spatial directions, and a time step of r = 5-10“"^. 

8.2. Numerical experiments. 

8.2.1. Qualitative behavior. In the first series of experiments, we choose the potential V as 

V{x) = -A( cos(27rxi) - l)( cos(47rx2) - l), (50) 

with A = 0.75. We approximate solutions for the initial density 

n° = C ( 0.1 + xi( cos(47rxi) - 1.2)( cos(27rj/2) “ l)) i 
where C > 0 is such that has unit mass. 



( 51 ) 
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Figure 1. Snapshots of the evolution of the fully discrete solution for the initial 
density (51), taken at times t = 0,2.5 • 10“^,4 • 10“^ and t = 5 - 10“^. Upper row: 
transport maps T^; Lower row: densities 


We employ our scheme with a relatively high spatial resolutions of K = 32. We visualize the 
result in Figure]^ the left column displays the transport maps T^ij the right column shows the 
corresponding density functions after various numbers n of time iterations. Qualitatively, the 
density evolves from the initial state that has two peaks of different height to the an approximation 
of the stationary solution that has two peaks of equal height, which are oriented in the orthogonal 
direction. Figure [^middle shows the correponding decay of the entropy E. This curve is compared 
to the one obtained for E(uref), using the reference solution with the same initial datum vP. The 
curves are virtually indistinguishable. 


8.2.2. Rate of convergence. For a quantitative estimate on the accuracy of our numerical scheme, 
we eliminate the external potential, F = 0, and run our scheme with the same initial datum 
as above in (51), with various different spatial resolutions, ranging from £ = 8 (coarsest) to £ = 24 
(finest). The time step width r = 5 • is kept fixed. 

These solutions are then compared to the reference solution Upef that has been obtained by the 
finite element scheme (49). In Figure [^left, we have plotted the behavior of the L^-norm of the 
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number K of modes 


time t 


time t ~ nr 


Figure 2. Left: L^-distance between the solution of our variational scheme and 
the reference solution at the final time T = 0.01, using i^T = 4,8,12,16,20,25 modes. 
Middle: Decay of the entropy. Right: Exponential rate of decay in time of the L^- 
distance between two particular fully discrete solutions. 


difference, evaluated at the time T = 0.01. The observed order of convergence is between two and 
three, i.e., the approximation error is roughly inverse proportional to the number of modes. 


8.2.3. Contractivity. Gradient flows of A-convex functionals are A-contractive. In the context at 
hand, this means that if the potential F in Q is A-convex, and ui, U 2 are two solutions to the 
initial value problem 0-(§ with respective initial data uf, U 2 , then 

m{ui{t),U2{t)) < W2{u°,4)e-^K (52) 


Since our discretization preserves A-convexity as detailed in Propostion a corresponding contrac¬ 
tion estimate holds — in each time step, see Corollary However, since the L^-space on which 
the minimization is performed changes from one time step to the next, it is not clear how to iterate 
that estimate and obtain an analogue of (52) above. 

A natural question to ask is if maybe a stronger contraction principle holds, namely if the L^- 
distance (with respect to u) between the transport maps for two solutions satisfies the analogous 
estimate as (521 for the Wasserstein distance between densities. In general, such a stronger principle 
is too much to ask for (recall that E is just poly-convex). However, some numerical experiments 
suggest that this estimate does hold if the initial data are already close in a suitable sense. 

For illustration, we have performed a numerical experiment for V{x) = ^||x ||2 with A = 10, with 
two initial conditions u? ^'^^^1 u® that are random perturbations of the unit density u = 1. More 
precisely, = (T°)#u for j = 1,2, with T?,*!® ^ ^<5 being random perturbations of the identity. 
The discretization parameters are AT = 12 and r = 10"^. The result is given in Figure [fright; there 
is apparently an exponential contraction of the transport maps, with an exponential rate that is 
higher (but still comparable to) A. 


8.2.4. Code. In principle, an actual implementation of our scheme can be quite compact as the 
conceptual code for the Julia language in Fig. shows. Note, however, that here we trade elegance 
for speed and that for the experiments reported on above we have actually used a much faster 
implementation. 
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using ForwardDiff 

d = 2; K = 5 

g(k,x) = cos(pi*k*x) 

dg(k,x) = -pi*k*sin(pi*k*x) 

ddg(k,x) = -pi'‘2*k“2* cos (pi*k*x) 

gradbCk,x) = 2/( pi *nor)n (k ) ) * [ dg ( k [ 1] , x [ 1] ) ♦ g (k [2] , x [2] ) ; gCk[l] 
Hb(k,x) = 2/Cpi*normCk))*[ddg(k [1] ,X [1])♦g(k [2] ,X [2]) dg(k[l],x 

dg(k[l] ,x[l])*dgCk[2] ,x[2]) gCk[l] 
I = [[kl,k2] for kl=0:K, k2=0:K] [2:end] 

tCz,x) = X + sum([z[k [2]*(K+1)+k [ 1]]*gradb(k,X) for k = I]) 
Dt(z,x) = eyeC2) + sum( [z [k [2]*(K+ 1)+k [ 1]]♦HbCk,x) for k = I]) 
hs(r) = 1/r 

L = 20; dx = 1/(L-1); dy = dx; 

X = [[x,y] for X = 0:dx:1, y = 0:dy:l] 
oml = ones(L); oml [1] = oml [L] = 1/2; 

om = [ a*b*dx''2 for a = oml , b = oml] ; 

u0(x) = (0.1 + x[l] * Ceos (4*pi*x[l])-l . 2) * Ceos C 2* pi *x [2] ) - 1) ) 

C = sumC[uOCx)*dx~2 for x = X]) 
ub = [uoex [k,1])/C for k = l:L, 1 = 1:L] 

VCx) = -0.75*CcosC2*pi*x[l])-l)*CcosC4*pi*x[2])-l) 
si = ones CL,L) 
tau = 5e-4 

intCf) = sumCf.*ub.♦om) 

WCz) = int C [sum C Ct Cz , X [k , 1] )-X [k , 1] ) .‘‘2) for k = l:L, 1 = 1:L]) 

ECz) = intC[ChsCsi[k,1]/ub[k,1]*detCDtCz,X[k,1]))) + VCtCz,X[k, 
FCz) = l/C2*tau)*WCz) + ECz) 
gradF = ForwardDiff.gradientCF) 

HF = ForwardDiff.hessianCF) 
z = zeros CsizeCI)); 
time = 0 

while time < 5*le-2 

while normCgradFCz)) > le-3 
dz = -HFCz)\gradFCz) 
z = z + dz 

end 

si = [det CDt Cz , X [k , 1] ) ) * si [k , 1] for k = 1:L, 1 = 1:L] 

X = [tCz,X[k,l]) for k = 1:L, 1 = 1:L]; 
time = time + tau 

end 


# dimension, number of modes 

,x[l])*dgCk[2] ,x[2])] 
[l])*dgCk[2] ,x[2]); 

. X [1] ) *ddg (k [2] , X [2] ) ] ; 

# Lagrangian maps 

# pressure function 

# space grid 

# weights for quadrature 


# reference density 

# potential 

# determinants 

# time step size 


1]))) for k = l:L , 1 = 1:L]) 


# time stepping 

# Newton iteration 


# updating determinants 

# updating space grid 


Figure 3. Conceptual code of the Lagrangian scheme in the Julia language. 


Appendix A. Gradient flows 

Lemma 18. Let two functionals E : 'P(n) ^ M and E : Diff^ ^ K &e given, which satisfy E(T) = 
E(ir^tt), for some reference density u. Then the derivatives of these functionals at any given 
T 6 Diff'*’, are related as follows: 

VJL 2 (d,) [DE(T#i2)] o T = [DE(T)]. (53) 

Proof. Let C ^ be a smooth vector field. Given T, consider the family (T®)s 6R of 

diffeomorphisms T® of Lt that is defined as solution to the ODE problem 

= T° = 5:. 

Be definition of E via E(T) = E(T^tt), we trivially have that 

E(T^u). (54) 

s=0 


_d 

ds 


E(T*) = 

T ds 
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We shall now compute both derivatives separately. On the one hand, 
d 


ds 


E(r) = DE(T) [C ° 2] = ^ {Jlhu) [DE(T)]} (0 ■ [C « T](e) duiO ■ 

s=0 


And on the other hand, with u 


ds 


E(T^u) = DE{u) 


s=0 


d 

ds 


s=0 


d 
’ ds 


dx 


= ;rl /'{JL^(d.)[DE(ii)]}or(0di2(0 

ls=0 

= ^ {vJi2(d,) [DE(ii)]} O TiO • [C o ^](0 di2(0- 


Recallin (541, and since ( was arbitrary, the equality ( |53[ ) follows. 

Lemma 19. The L‘^{%^u)-gradient of the functional E defined in (16) is given by 


□ 


v(t;T)=;^div 


IdetDtj ^ ^ 


+ VEot. 


Proof. Fix t 6 VX^. From the definition in (16), we obtain for every “variation” ^ € VTX and 
sufficiently small e: 

” b ■ /„ {'*. " 'o} 

Now take the derivative with respect to £ at £ = 0. For simplification of the result, use that 
h'^{s) = -P(l/s), and perform an integration by parts: 


d£ 


* ^0 “d ■ /„ (1^) td(D*)-‘Dp. f. vrd)} iT,f. 




dy + y~ ‘C'VE(i) d%^u 

+ S/V o • C dT^u. 


This yields the result. 


□ 
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